/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */
#include "getformatXML.h"

#include "l1_olci.h"
#include "l1c_msi.h"
#include <genutils.h>
#include <pugixml.hpp>
#include <stdlib.h>
#include <string.h>
#include <timeutils.h>
#include <boost/foreach.hpp>
#include <vector>
#include <cmath>
#include <cfloat>

#include "l1c_msi_private.h"
#include <boost/algorithm/string.hpp>


using namespace pugi;
using namespace std;

/**
 *   chk_olci
 *   purpose: check a file to see if it is an OLCI xml file
 *   If it is, populate the data pointer with filenames
 *   associated with Lt's, instrument and geometry data
 *   Returns -1 if not OLCI or the format code
 *
 * @param char *fname      - name of file to check
 * @param filehandle *file - input file information
 *
 * Author: R. Healy (SAIC) 7/29/2015
 * W. Robinson, SAIC, 18 Feb 2016 pad out the olci file name storage by 1 
 *   to accomodate trailing null, allow xml files with alternate product 
 *   name to be accepted as olci files
 * s
 */

extern "C" int32_t chk_olci_xml(char *fname, filehandle *file) {
    xml_document rootNode;
    xml_node dataNode;
    olci_t *data = NULL;

    const char *productName, *dataName;
    int nbands = 0;
    int32_t fmt = -1;

    if (!rootNode.load_file(fname)) {
        return fmt;
    }

    productName = rootNode.first_element_by_path("xfdu:XFDU/informationPackageMap/xfdu:contentUnit").attribute("textInfo").value();
    if (strstr(productName, "OLCI") &&
            strstr(productName, "SENTINEL-3") &&
            strstr(productName, "Level 1")) {
        if (want_verbose)
            printf("%s\n", productName);
        file->private_data = data = createPrivateData_olci();
        fmt = FMT_OLCI;
    } else {
        return fmt;
    }

    int foundGeoCoord = 0;
    int foundtieGeoCoord = 0;
    int foundtieGeometries = 0;
    int foundinstrumentDat = 0;
    int foundtimeCoord = 0;
    int foundTieMeteo = 0;
    
    string indir = file->name;
    size_t index = indir.find_last_of('/');
    if(index != string::npos)
        indir.erase(index);
    else
        indir.clear();
    string fileName;
    
    dataNode = rootNode.first_element_by_path("xfdu:XFDU/dataObjectSection").first_child();
    while (dataNode) {
        dataName = dataNode.attribute("ID").value();
        productName = dataNode.first_element_by_path("byteStream/fileLocation").attribute("href").value();
        if(indir.empty())
            fileName = productName;
        else
            fileName = indir + "/" + productName;
        index = fileName.find("./");
        if(index != string::npos)
            fileName.erase(index, 2);
        
        if (strstr(dataName, "radianceData")) {
            if (nbands > MAXOLCI_RADFILES) {
                printf("%s, %d - E - Maximum number of radiance files (%d) reached\n",
                        __FILE__, __LINE__, MAXOLCI_RADFILES);
                exit(EXIT_FAILURE);
            }
            data->olci_radfiles[nbands] = strdup(fileName.c_str());
            if (!data->olci_radfiles[nbands]) {
                printf("%s, %d - E - unable to allocate radiance filename \n",
                        __FILE__, __LINE__);
                exit(EXIT_FAILURE);
            }
            nbands++;
        } else if (strcmp(dataName, "geoCoordinatesData") == 0) {
            data->geoCoordinatesFile = strdup(fileName.c_str());
            if (!data->geoCoordinatesFile) {
                printf("%s, %d - E - unable to allocate geo coordinates filename \n",
                        __FILE__, __LINE__);
                exit(EXIT_FAILURE);
            }
            foundGeoCoord = 1;
        } else if (strcmp(dataName, "tieGeoCoordinatesData") == 0) {
            data->tieGeoCoordinatesFile = strdup(fileName.c_str());
            if (!data->tieGeoCoordinatesFile) {
                printf("%s, %d - E - unable to allocate tie geo coordinates filename \n",
                        __FILE__, __LINE__);
                exit(EXIT_FAILURE);
            }
            foundtieGeoCoord = 1;
        } else if (strstr(dataName, "tieGeometriesData")) {
            data->tieGeometriesFile = strdup(fileName.c_str());
            if (!data->tieGeometriesFile) {
                printf("%s, %d - E - unable to allocate tie geometries filename \n",
                        __FILE__, __LINE__);
                exit(EXIT_FAILURE);
            }
            foundtieGeometries = 1;
        } else if (strstr(dataName, "instrumentDataData")) {
            data->instrumentFile = strdup(fileName.c_str());
            if (!data->instrumentFile) {
                printf("%s, %d - E - unable to allocate instrument filename \n",
                        __FILE__, __LINE__);
                exit(EXIT_FAILURE);
            }
            foundinstrumentDat = 1;
        } else if (strstr(dataName, "timeCoordinatesData")) {
            data->time_coordinatesFile = strdup(fileName.c_str());
            if (!data->time_coordinatesFile) {
                printf("%s, %d - E - unable to allocate time_coordinates filename \n",
                        __FILE__, __LINE__);
                exit(EXIT_FAILURE);
            }
            foundtimeCoord = 1;
        } else if (strstr(dataName, "tieMeteoData")) {
            data->tieMeteoFile = strdup(fileName.c_str());
            if (!data->tieMeteoFile) {
                printf("%s, %d - E - unable to allocate tieMeteoFile filename \n",
                        __FILE__, __LINE__);
                exit(EXIT_FAILURE);
            }
            foundTieMeteo = 1;
        }
        dataNode = dataNode.next_sibling();
    }
    if (foundtimeCoord && foundinstrumentDat && foundtieGeometries && foundtieGeoCoord && foundtieGeoCoord && foundGeoCoord && foundTieMeteo && nbands == MAXOLCI_RADFILES) {
        if (want_verbose)
            printf("OLCI:  Found all meta data information\n");
    } else {
        printf("%s, %d - E - missing meta data information in file %s \n",
                __FILE__, __LINE__, fname);
        exit(EXIT_FAILURE);
    }
    return fmt;
}

extern "C" int32_t chk_msi_xml(char *filename, filehandle *file) {
    xml_document rootNode;
    xml_node dataNode, metaNode;
    msi_t *data = NULL;

    const char *productName, *dataName;
    int nbands=0;
    int32_t year, doy, msec;
    int32_t fmt = -1;

    if(!rootNode.load_file(filename)) {
       return fmt;
    }
    
    productName = rootNode.first_element_by_path("xfdu:XFDU/informationPackageMap/xfdu:contentUnit").attribute("textInfo").value();
    if ( strstr(productName,"MSI" ) &&
        strstr(productName,"SENTINEL-2" )  &&
        strstr(productName,"Level-1C" ) ) {
        if(want_verbose)
            printf("Found Sentinel-2 product: %s\n",productName);

        fmt = FMT_MSIL1C;
        file->format = fmt;
        file->sensorID = MSIS2A;
        file->subsensorID = MSI_S2A;
        // Allocate space for private data
        file->private_data = data = createPrivateData_msi(maxBands); 
    } else {
        return fmt;
    }
    
    metaNode = rootNode.first_element_by_path("xfdu:XFDU/metadataSection");

    // platform ID
    dataNode = metaNode.find_child_by_attribute("metadataObject", "ID", "platform");
    dataName = dataNode.first_element_by_path("metadataWrap/xmlData/safe:platform/safe:number").child_value();
    if(strcmp(dataName, "2B") == 0) {
        file->sensorID = MSIS2B;
        file->subsensorID = MSI_S2B;
    }
        
    // orbit direction
    dataNode = metaNode.find_child_by_attribute("metadataObject", "ID", "measurementOrbitReference");
    dataName = dataNode.first_element_by_path("metadataWrap/xmlData/safe:orbitReference/safe:orbitNumber").attribute("groundTrackDirection").value();
    data->orbit_direction = strdup(dataName);

    // Band image paths
    int nbandsImg = 0;
    int nbandsDetFoot = 0;
    int pickMe = 0;
    
    dataNode = rootNode.first_element_by_path("xfdu:XFDU/dataObjectSection/dataObject") ;

    string indir = file->name;
    size_t index = indir.find_last_of('/');
    if(index != string::npos)
        indir.erase(index);
    else
        indir.clear();
    string fileName;
    
    while (dataNode) {

        dataName = dataNode.attribute("ID").value();
        productName = dataNode.first_element_by_path("byteStream/fileLocation").attribute("href").value();
        if(indir.empty())
            fileName = productName;
        else
            fileName = indir + "/" + productName;
        index = fileName.find("./");
        if(index != string::npos)
            fileName.erase(index, 2);

        if (strstr(dataName, "S2_Level-1C_Tile1_Metadata")){
            data->granuleMetadataFile = strdup(fileName.c_str());
        }
        if (strstr(dataName, "S2_Level-1C_Datastrip1_Metadata")){
            data->datastripMetadataFile = strdup(fileName.c_str());
        }
        if (strstr(dataName, "IMG_DATA_Band") && !(strstr(dataName, "TCI"))) {
            if (nbandsImg > maxBands) {
                printf("%s, %d - E - Maximum number of radiance values (%d) reached\n",
                        __FILE__, __LINE__, maxBands);
                exit(EXIT_FAILURE);
            }
            data->jp2[nbandsImg] = strdup(fileName.c_str());
            if (!data->jp2[nbandsImg]) {
                printf("%s, %d - E - unable to set path for band %d\n",
                        __FILE__, __LINE__, nbandsImg);
                exit(EXIT_FAILURE);
            }
            nbandsImg++;
        }
        if (strstr(dataName, "DetectorFootprint")) {
            if (nbandsDetFoot > maxBands) {
                printf("%s, %d - E - Maximum number of radiance values (%d) reached\n",
                        __FILE__, __LINE__, maxBands);
                exit(EXIT_FAILURE);
            }
            data->detectorFootprintFiles[nbandsDetFoot] = strdup(fileName.c_str());
            if (!data->detectorFootprintFiles[nbandsDetFoot]) {
                printf("%s, %d - E - unable to set path for band %d\n",
                        __FILE__, __LINE__, nbandsDetFoot);
                exit(EXIT_FAILURE);
            }
            if (strstr(data->detectorFootprintFiles[nbandsDetFoot], "B07")) {
                pickMe = nbandsDetFoot;
            }
            nbandsDetFoot++;    
        }
        dataNode = dataNode.next_sibling("dataObject");
    }
    
    /* Read tile metadata file */
    if(!readTileMeta_msi(file))
        fprintf(stderr, "-E- %s line %d: unable read tile metadata file for MSI dataset %s\n",
            __FILE__,__LINE__,data->granuleMetadataFile);
    
    /* Read datastrip metadata file */
    if(!readDatastripMeta_msi(file))
        fprintf(stderr, "-E- %s line %d: unable read datastrip metadata file for MSI dataset %s\n",
            __FILE__,__LINE__,data->datastripMetadataFile);
    
    /* Read detector footprint file 
     for simplicity, read just one for a 20m band
     */
    if(!readDetectorFootprint_msi(file, pickMe))
        fprintf(stderr, "-E- %s line %d: unable read detector footprint file for MSI dataset %s\n",
            __FILE__,__LINE__,data->detectorFootprintFiles[pickMe]);

    return fmt;
}

int32_t readTileMeta_msi(filehandle *file) {
    xml_document rootNode;
    xml_node dataNode;
    xml_node zenithNode;
    xml_node azimuthNode;
    xml_node valuesNode;
    char *delim = " "; // input separated by spaces
    char *token = NULL;
    msi_t *data = (msi_t*) file->private_data;
    char* tmpBuff;

    const char *xmlBuffer, *dataName;
    if(!rootNode.load_file(data->granuleMetadataFile)) {
       return 0;
    }
    
    // fix orbit altitude
    data->alt = 786.;
    // Get geocoding information
    dataNode = rootNode.first_element_by_path("n1:Level-1C_Tile_ID/n1:Geometric_Info/Tile_Geocoding");

    // Get coordinate system EPSG code
    xmlBuffer = dataNode.first_element_by_path("HORIZONTAL_CS_CODE").first_child().value();
    std::string s = xmlBuffer;

    std::string delimiter = ":";
    s.erase(0, s.find(delimiter) + delimiter.length());
    data->CSCode = atoi(s.c_str());
    
    // Get coordinate system UTM zone

    xmlBuffer = dataNode.first_element_by_path("HORIZONTAL_CS_NAME").first_child().value();
    s = xmlBuffer;

    delimiter = "zone ";
    s.erase(0, s.find(delimiter) + delimiter.length());
    data->UTMZone = strdup(s.c_str());
       
    xmlBuffer = dataNode.find_child_by_attribute("Size","resolution","20").child_value("NCOLS");
    file->npix = atoi(xmlBuffer);
    xmlBuffer = dataNode.find_child_by_attribute("Size","resolution","20").child_value("NROWS");
    file->nscan = atoi(xmlBuffer);
    xmlBuffer = dataNode.find_child_by_attribute("Geoposition","resolution","20").child_value("ULX");
    data->ULCoord[0] = atoi(xmlBuffer);    
    xmlBuffer = dataNode.find_child_by_attribute("Geoposition","resolution","20").child_value("ULY");
    data->ULCoord[1] = atoi(xmlBuffer);    
    

    // Set projection string
    char pjStr[FILENAME_MAX];
    sprintf(pjStr, "+proj=utm +ellps=WGS84 +datum=WGS84 +zone=%s +units=m", data->UTMZone);

    if (!(data->pj_latlong = pj_init_plus("+proj=longlat +datum=WGS84"))) {
        fprintf(stderr, "-E- %s line %d: pj_init_plus error: longlat \n",
                __FILE__, __LINE__);
        exit(EXIT_FAILURE);
    }
    if (!(data->pj_utm = pj_init_plus(pjStr))) {
        fprintf(stderr, "-E- %s line %d: pj_init_plus error: utm \n",
                __FILE__, __LINE__);
        exit(EXIT_FAILURE);
    }
    
    // Read tie-point sensor view angles
    dataNode = rootNode.first_element_by_path("n1:Level-1C_Tile_ID/n1:Geometric_Info/Tile_Angles/Viewing_Incidence_Angles_Grids");

      // Set up tie point array sizes ([numDetectors][maxBands][HEIGHT][WIDTH])
    int tiepoint_height = 23;
    int tiepoint_width = 23;

    data->sensorZenith = (double **) malloc(maxBands * sizeof (double *));
    data->sensorAzimuth = (double **) malloc(maxBands * sizeof (double *));
    for (int i = 0; i < maxBands; i++) {
        data->sensorZenith[i] = (double *) calloc(tiepoint_width * tiepoint_height, sizeof (double ));
        data->sensorAzimuth[i] = (double *) calloc(tiepoint_width * tiepoint_height, sizeof (double ));
    }

    while (dataNode) {
        xmlBuffer = dataNode.attribute("bandId").value();
        int bandIdx = atoi(xmlBuffer);
        xmlBuffer = dataNode.attribute("detectorId").value();
        int detectorIdx = atoi(xmlBuffer) - 1;
        zenithNode = dataNode.first_element_by_path("Zenith/Values_List");
        valuesNode = zenithNode.first_element_by_path("VALUES");
        int i = 0;
        while (valuesNode){
            tmpBuff = strdup(valuesNode.first_child().value());
            int j = 0;
            for (token = strtok((char *) tmpBuff, delim); token != NULL; token = strtok(NULL, delim)) {
                char *unconverted;
                double value = strtof(token, &unconverted);
                if (!std::isnan(value))
                    data->sensorZenith[bandIdx][i*tiepoint_width + j] = value;
                j++;
            }
            free(tmpBuff);
            i++;
            valuesNode = valuesNode.next_sibling("VALUES");
        }
        azimuthNode = dataNode.first_element_by_path("Azimuth/Values_List");
        valuesNode = azimuthNode.first_element_by_path("VALUES");
        i = 0;
        while (valuesNode){
            tmpBuff = strdup(valuesNode.first_child().value());
            int j = 0;
            for (token = strtok((char *) tmpBuff, delim); token != NULL; token = strtok(NULL, delim)) {
                char *unconverted;
                double value = strtod(token, &unconverted);
                if (!std::isnan(value))
                    data->sensorAzimuth[bandIdx][i*tiepoint_width + j] = value;
                j++;
            }
            free(tmpBuff);
            i++;
            valuesNode = valuesNode.next_sibling("VALUES");
        }
        dataNode = dataNode.next_sibling("Viewing_Incidence_Angles_Grids");
    }
    return 1;
}

int32_t readDatastripMeta_msi(filehandle *file) {
    xml_document rootNode;
    xml_node dataNode;
    xml_node subDataNode;
    xml_node_iterator it;
    double unixseconds;
    char *delim = " "; // input separated by spaces
    char *token = NULL;
    msi_t *data = (msi_t*) file->private_data;
    char* tmpBuff;

    const char *xmlBuffer, *dataName;
    if (!rootNode.load_file(data->datastripMetadataFile)) {
        return 0;
    }

    // Get start time
    dataNode = rootNode.first_element_by_path("n1:Level-1C_DataStrip_ID/n1:General_Info/Datastrip_Time_Info");

    xmlBuffer = rootNode.first_element_by_path("n1:Level-1C_DataStrip_ID/n1:General_Info/Datastrip_Time_Info/DATASTRIP_SENSING_START").first_child().value();
    data->scene_start_time = isodate2unix((char*) xmlBuffer);

    // Get detector line period
    dataNode = rootNode.first_element_by_path("n1:Level-1C_DataStrip_ID/n1:Image_Data_Info/Sensor_Configuration/Time_Stamp");
    xmlBuffer = dataNode.first_element_by_path("LINE_PERIOD").first_child().value();
    data->lineTimeDelta = atof(xmlBuffer) * 2e-3; //convert to seconds, for 20m band
    
    // Get detector relative startimes
    dataNode = rootNode.first_element_by_path("n1:Level-1C_DataStrip_ID/n1:Image_Data_Info/Sensor_Configuration/Time_Stamp/Band_Time_Stamp");
    while (dataNode) {
        xmlBuffer = dataNode.attribute("bandId").value();
        int bandIdx = atoi(xmlBuffer);
        subDataNode = dataNode.first_element_by_path("Detector");
        while (subDataNode) {
            xmlBuffer = subDataNode.attribute("detectorId").value();
            int detectorIdx = atoi(xmlBuffer) - 1;
            xmlBuffer = subDataNode.child_value("GPS_TIME");
            unixseconds = isodate2unix((char*) xmlBuffer);
            data->detectorDeltaTime[detectorIdx][bandIdx] = data->scene_start_time - unixseconds;
            subDataNode = subDataNode.next_sibling("Detector");
        }
        
        dataNode = dataNode.next_sibling("Band_Time_Stamp");
    }
    
    // Get number of GPS entries
    dataNode = rootNode.first_element_by_path("n1:Level-1C_DataStrip_ID/n1:Satellite_Ancillary_Data_Info/Ephemeris/GPS_Points_List");

    data->num_gps = 0;
    for (it = dataNode.begin(); it != dataNode.end(); it++) {
        data->num_gps++;
    }
    data->gpstime = (double *) malloc(data->num_gps * sizeof(double));
    for (int i = 0; i<3; i++)
        data->position[i] = (double *) malloc(data->num_gps * sizeof(double));

    // read GPS
    int i = 0;

    dataNode = rootNode.first_element_by_path("n1:Level-1C_DataStrip_ID/n1:Satellite_Ancillary_Data_Info/Ephemeris/GPS_Points_List/GPS_Point");

    while (dataNode) {
        tmpBuff = strdup(dataNode.child_value("POSITION_VALUES"));
        int j = 0;
        for (token = strtok((char *)tmpBuff, delim); token != NULL; token = strtok(NULL, delim)) {
            char *unconverted;
            data->position[j][i] = strtod(token, &unconverted);
            j++;
        }
        free(tmpBuff);
        xmlBuffer = dataNode.child_value("GPS_TIME");
        data->gpstime[i] = isodate2unix((char*) xmlBuffer);

        i++;
        dataNode = dataNode.next_sibling("GPS_Point");
    }

    return 1;
}

int32_t readDetectorFootprint_msi(filehandle *file, int band) {

    xml_document rootNode;
    xml_node dataNode;
    xml_node polyNode;
    msi_t *data = (msi_t*) file->private_data;
    std::vector<std::string> detstrs;
    std::vector<std::string> polypointstrs;
    std::vector<std::string>::iterator sit;

    const char *xmlBuffer, *detectorName, *polygon;
    if (!rootNode.load_file(data->detectorFootprintFiles[band])) {
        return 0;
    }

    dataNode = rootNode.first_element_by_path("eop:Mask/eop:maskMembers/eop:MaskFeature");
    while (dataNode) {
        detectorName = dataNode.attribute("gml:id").value();
        //detector_footprint-B05-03-0
        boost::split(detstrs, detectorName, boost::is_any_of("-"));
        int detidx = std::atoi(detstrs[2].c_str()) - 1;

        polyNode = dataNode.first_element_by_path("eop:extentOf/gml:Polygon/gml:exterior/gml:LinearRing/gml:posList");
        polygon = polyNode.first_child().value();

        boost::split(polypointstrs, polygon, boost::is_any_of(" "));
        std::string polyWKT = "POLYGON((";
        int i = 1;
        for(sit=polypointstrs.begin() ; sit < polypointstrs.end(); sit++,i++ ) {
        // skip every third element
            if(i % 3) 
                polyWKT = polyWKT + *sit + " ";
            else
                polyWKT = polyWKT + ",";

        }
        polyWKT = polyWKT + "))";
        boost::geometry::read_wkt(polyWKT, data->detectorFootprints[detidx]);

        dataNode = dataNode.next_sibling("eop:MaskFeature");
    }
    
    return 1;
}
